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Abstract 

We focus on the spatial discretization produced by the Variational Particle-Mesh (VPM) 
method for a prototype fluid equation the known as the EPDiff equation, which is short 
for Euler-Poincare equation associated with the diffeomorphism group {of R'*, or of a d- 
dimensional manifold fi). The EPDiff equation admits measure valued solutions, whose 
dynamics are determined by the momentum maps for the left and right actions of the dif- 
feomorphisms on embedded subspaces of R''. The discrete VPM analogs of those dynamics 
are studied here. Our main results are: (i) a variational formulation for the VPM method, 
expressed in terms of a constrained variational principle principle for the Lagrangian particles, 
whose velocities are restricted to a distribution 7?vpm which is a finite-dimensional subspace 
of the Lie algebra of vector fields on 51; {ii) a corresponding constrained variational principle 
on the fixed Eulerian grid which gives a discrete version of the Euler-Poincare equation; and 
{iii) discrete versions of the momentum maps for the left and right actions of diffeomorphisms 
on the space of solutions. 
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1 Introduction 

1.1 Transverse internal wave interactions 

Synthetic Aperture Radar (SAR) observations from the Space Shuttle often show nonlinear internal 
wave trains that propagate for many hundreds of kilometers across large basins such as the South 
China Sea (SCS) shown in Figure Q 

These wave trains are characterized as Great Lines on the Sea in |16| . Both lines and spirals 
on the sea arise as flow phenomena, rather than wave phenomena per se. The flow phenomenon 
detected in the the SAR imagery is associated with nonlinear internal waves, whose crests may be 
as much as 200km long. The amplitude of these internal waves results in about 150m of deflection 
in the thermocline over a distance of about 1 km. Thus, their aspect ratio satisfies the first criterion 
to be nonlinear shallow water waves. Their amplitude is also considerably less than the typical 
thickness of the thermocline, but it is not actually infinitesimal compared to the thermocline 
thickness. The flow along the crests of these waves also indicates they are not precisely the same 
as usual shallow water waves. 

The particular nonlinear internal waves found in the SCS are generated by the tides flowing 
East to West through the Luzon Strait over submerged ridges between Taiwan and the Phillipines. 
The SAR images in Figure ^ show that the momentum of the tides flowing Westward over these 
ridges concentrates into internal waves on the thermocline that emerge into the SCS basin as thin 
wave fronts which may extend in length for hundreds of kilometers (much larger than the Straits in 
which they were created) and may propagate for thousands of kilometers. Perhaps because of the 
complex topography, the tides flowing over the mouth of the Luzon Strait do not produce internal 
waves propagating in both directions. The signiflcant wave trains propagate Westward. 

Propagating wave trains may intersect transversely with other wave trains. Sometimes these 
wave trains merely pass through each other as linear waves. However, in nonlinear wave encounters 
such as those captured by SAR imaging of the region of the SCS West of Dong Sha Island in Figure 
121 two wave fronts may intersect transversely, merge together and produce a single wave front. This 
merger of the wave fronts is the hallmark of a nonlinear process. These particular wave interactions 
possess strong transverse dynamics (flow along the crests) and momentum exchange in the direction 
of propagation, which allow the wave fronts to merge and reconnect, rather than merely passing 
through each other, as weaker waves do when they intersect in an interference pattern. 

Nonlinear internal wave interactions have been well studied in one dimension, often by using 
the weakly nonlinear Boussinesq approximation. These studies have usually resulted in a variant 
of the Korteweg-de Vries (KdV) equation, which has soliton solutions that interact by exchange 
of momentum in unidirectional elastic collisions (Whitham 1967). However, the complex wave 
front interactions shown in in Figure [3 are plainly at least two-dimensional. We shall pursue 
the qualitative description of these higher-dimensional wave interactions by using a simple two- 
dimensional model equation called EPDiff.^ EPDiff may be derived in one dimension from the 
asymptotic expansion for shallow water wave motion of the Euler equations for the unidirectional 
flow of an incompressible fluid with a free surface moving under gravity. In one dimension the result 
is the Camassa-Holm (CH) equation, which arises at quadratic order in this expansion. That is, 
CH is one order of accuracy in the asymptotic expansion beyond KdV, which arises at linear order. 
Just as for KdV, the CH equation is completely integrable; so CH also has soliton solutions that 
interact by elastic collisions in one dimension. Moreover, in the limit of zero linear dispersion, the 
CH solitons develop a sharp peak at which their profile has a jump in derivative that forms a sharp 
peak. In this limit, the CH solitons are called "peakons." The CH peakons are weak solutions, in 
the sense that their momentum is concentrated on delta functions that move with the velocity of 
the fluid flow. 

^EPDiff is the "Euler- Poincarc equation on the diflfeomorphisms" . 
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In its zero-dispersion limit, CH has a geometric property that aUows it to be immediately gen- 
eralized to higher dimensions, in which it is called EPDifF. The term "EPDifF" distinguishes CH, 
which is a one-dimensional shallow water wave equation with physical wave dispersion, from its 
dispersionless limit which belongs to a larger class of equations. This larger class of equations - the 
Euler-Poincare (EP) equations ~- describes geodesic motion with respect to any metric defining 
a norm on the vector space of the Lie algebra of a Lie group. In the geometric theory of fluid 
mechanics, the fluid velocity belongs to the tangent space of the group of smooth invertible maps, 
called "diffcomorphisms" (or diffcos, for short). The Euler-Poincare equation on the diffcomor- 
phisms is called EPDiff. EPDiff is a larger class of equations than CH also because it is defined for 
geodesic motion on the diffeos with respect to any metric, not just for the norm of the velocity, 
which appears as the kinetic energy norm in the derivation of CH. (The gradient part of the 
norm for CH corresponds to the vertically averaged kinetic energy associated with vertical motion.) 
Thus, among the EPDiff equations, the dispersionless limit of CH is one-dimensional EPDiff(i7^). 
In one dimension, the momentum of the EPDiff(i7^) peakons is concentrated at points moving 
along with the flow; but in higher dimensions, their momentum is distributed on embedded sub- 
spaces moving with the flow. In particular, EPDiff (i7^) in two dimensions has singular solutions 
whose momentum is distributed along curves in the plane. As solutions of the two-dimensional 
version of a unidirectional shallow water wave equation in its limit of zero linear dispersion, these 
moving curves in the plane evolving under the dynamics of EPDiff(i/^) arc prototypes for studying 
the interactions of the Great Lines on the Sea. 

To jump ahead, the singular (or, weak) solutions of the EPDiff equation that emerge in finite 
time from any confined smooth initial conditions and are supported on embedded subspaces moving 
with the flow velocity, just as seen in the Great Lines on the Sea captured in Figure We 
developed a numerical method for simulating the singular solutions of EPDiff in the framework of 
its geometric definition, which is natural for the Variational Particle Mesh (VPM) method. Our 
numerical results using VPM show that 

• Singular solutions for EPDiff may be simulated by VPM as curve-segments moving with the 
2D flow velocity that possess no internal degrees of freedom. 

• In collisions between any two of these curve-segment solutions for EPDiff, the momentum of 
the one that overtakes from behind is imparted to the one ahead. Thus, overtaking collisions 
between two finite-length wave packets are elastic. 

• The transverse collision of two curve-segment solutions for EPDiff may result in merger (or, 
reconnection) of the curve segments due to a combination of exchange of momentum between 
the wave trains and ffow along their wave crests. In two dimensions, the reconnection or 
merger of singular wave fronts under numerical EPDiff dynamics using VPM is evident in 
Figure El 

Plan of the paper In this paper we introduce the VPM method for EPDiff, and discuss some 
of the properties that arise from the variational structure, in the following sections: 

• The particle- mesh calculus is set out in section |5] 

• We give a variational principle associated with the method in section |3| 

• Section ^ shows that the Eulerian grid quantities satisfy an approximation the the EPDiff 
equation in Euler-Poincare form. 

• Section El defines a left action of I^vpm on fl and a provides the corresponding momentum 
map. 

• Sections El defines a right action in an extended space which can be interpreted as a discrete 
form of relabelling of Lagrangian particles. The Hamiltonian for the continuous time evo- 
lution of discretiscd EPDiff solutions is invariant under the action and so from Noether's 
theorem wc obtain a conserved momentum. 
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• Section [T] shows how this conserved momentum can be interpreted as a discrete form of 
Kelvin's circulation theorem. 

• Section |S1 gives some numerical examples, as well as convergence tests for the method. 




Figure 1: Image from the Space Shuttle of long, tidally-excited waves in the South China Sea near 
the Dong sha Atoll. The waves are propagating from East to West, and are produced every tide 
(about 12 hours). The waves interact with the Atoll and then undergo nonlinear reconnections. 
Picture courtesy of A. Liu. 



1.2 Theoretical development 

Much of the theoretical development in this paper is inspired by the following theorem 

Theorem 1.1 (Arnold (1966) U)- The solutions of Euler's equations for the incompressible 
motion of an ideal fluid describe coadjoint geodesic motion on the volume preserving diffeomor- 
phisms, with respect to the norm of the fluid velocity (the kinetic energy). 

The Euler equations for incompressible motion of an ideal fluid may be written in the material 
frame as 

„ / du \ , dx , „ 

"P — =0 along — = u with V ■ u = , 
\ dt J dt 

where V is the Leray projection onto the incompressible vector fields. These equations may also 
be written in the spatial frame as 

■p(<9tu + ad>) =0, 

where ad* is the dual of the ad-action among incompressible vector fields under the pairing. 
That is, ad*is defined by (/x , aduw) = — (ad*/x, w). Here ad„«; = [u,w] is the Lie-algebra 
commutator between vector fields u, w, and ( • , • ) denotes the pairing between such vector 
fields and one-form densities such as fi. 

EPDiff 
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Figure 2: Enlargement of part of figure ^ showing reconnecting long waves. Picture courtesy of 
A. Liu. 

The EPDiff equation describes the corresponding coadjoint geodesic motion on the full diffeo- 
morphism group, allowing for compressibility and an arbitrary norm. || • ||, 

S£ 1 
EPDiff is dtH + a.d*fi = , with fi = — , where £ = -\\uP . 

ou 2 

The momentum density /x is a one-form density and the EPDiff equation describes coadjoint 
dynamics under the action of the corresponding velocity vector field. In EPDiff, ad^ is the coadjoint 
action of a vector field u acting on a one-form density ^ 5£/Su for a Lagrangian £[u\ in Hamilton's 
principle dS = for S = J £[u\dt. In components, 

/X = m • dx (g) dVol 

and EPDiff may be written as the invariancc condition, 

—— = along — = u = G * m , 
at at 

where G* denotes convolution with the Green's function relating the components of m and u. In 
particular, for the norm ||ttp = / |up + q;^|Vu|^ dVol, we have the component relation 

m = u-Q;^Au, (1) 

and G is the Green's function for the Helmholtz operator. Id — A, A is the Laplacian, and a is 
a lengthscale. Thus, EPDiff for the H'^ norm with s > is an integro-partial differential equation. 

Originally derived as an n-dimensional generalisation of the Camassa-Holm equation for 
shallow- water dynamics in one dimension [2], EPDiff arises in several other applications. For 
example, EPDiff for the norm is the pressureless version of the Lagrangian-averaged Navier- 
Stokes-alpha (LANS-alpha) model of turbulence jH]. EPDiff for also emerges in the limit in 
which one ignores variations in height of the Green-Nagdhi equation for shallow water dynamics 
[3]. In one dimension, this is the dispersionless limit of the Camassa-Holm equation |2- In general, 
EPDiff is the equation for coadjoint geodesic motion on the diffeomorphisms with repect to any 
given norm on the Eulerian particle velocity (kinetic energy). Finally, EPDiff also describes the 
process of template matching in computational anatomy |15j . In this application, EPDiff has 
recently become a conduit for technology transfer from soliton theory to computational anatomy 
|12| . Thus, EPDiff turns out to be a prototype equation for a number of applications. 

The present article describes the underlying principles for using the Variational Particle-Mesh 
(VPM) method in numerically integrating EPDiff in the study of its nonlinear wave interactions. 
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Variational Particle-Mesh (VPM) method The Variational Particle-Mesh (VPM) method 
introduced in |S] produces Hamiltonian spatial discretizations of fluid equations which may then be 
integrated in discrete time by using a variational integrator. VPM may be regarded as a descendant 
of the Hamiltonian Particle-Mesh method , which is a Hamiltonian discretization of the rotating 
shallow-water equations. The difference is that HPM combines an Eulerian representation of the 
potential energy (which gives rise to the pressure term) with a Lagrangian representation of the 
kinetic energy, whilst VPM uses an Eulerian representation of the entire Lagrangian. This means 
that the VPM method is much more general than HPM and may be applied to many different fluid 
PDEs {e.g., shallow- water, Green-Nagdhi, incompressible Euler, etc.). In this paper we focus on 
EPDiff, which is an equation for fluid velocity only. Consequently, symmetries of the discretised 
fluid velocity will be symmetries of the equations. In future we will extend this work to include 
advected quantities such as density, scalars etc. Our ultimate aim is to use geometric properties 
in constructing general numerical methods for PDEs describing the continuum dynamics of fluids, 
complex fluids and plasmas. 

The conservative properties of variational integrators are well understood |18| . In this article, 
we will discuss preservation under VPM spatial discretization of the geometric properties of the 
well-known EPDiff equation for coadjoint motion under the diffeomorphisms |1(J| . 

9t/i + adlfi = , with fi = 

ou 

In particular, we shall discuss discrete VPM analogs of the momentum maps for the left and right 
actions of the diffeomorphisms on embedded subspaces of JOI ■ The Lagrangian we shall choose 
is the norm, i[u] = ^ ||i6||^i , so the components of velocity u and momentum density /x = blfbu 
will be related by the Helmholtz operator, as in equation In this case, velocity u S iJ^ implies 
that its dual momentum density fi G so the solutions of EPDiff may be measure valued in 

/x. That is, weak solutions of EPDiff are allowed in this case, which are expressed in terms of delta 
functions supported on embedded subspaces of R" ^Hl • The left action of the diffeomorphisms on 
these embedded subspaces of R" generates the motion of spatially discrete EPDiff (lattice EPDiff), 
while the right action is a symmetry and generates the conservation law for circulation according to 
the Kelvin-Noether theorem Thus, we seek the spatially discrete version of the corresponding 
theorem for continuum solutions in ^U). All of these properties will then be preserved by an 
appropriate variational time integrator. 

2 Particle-mesh calculus 

This section describes the particle-mesh calculus that will be used in discretising EPDiff. We 
shall describe its discretisation in space with continuous time, and later we shall describe how to 
construct variational time integrators to assemble a fully discrete space-time integration scheme. 

A finite dimensional subspace of X(il) The infinite-dimensional space of smooth vector fields 
X(f2) generates the diffeomorphisms (smooth invcrtible maps with smooth inverse) of the domain 
f2 onto intself. To make a numerical algorithm that can be calculated on a finite computer, we first 
need to choose a finite-dimensional subspace Xo of X(0) that will generate our diffeomorphisms. We 
begin with a fixed grid consisting of Ug points in the domain f2 with vector coordinates {fJfelfc^,! G 
j^dxrig ^ dimensions. At each grid point we shall associate a velocity vector u^. The 
finite-dimensional space of possible sets of velocity vectors u = (mi, . . . ,M,ig) G R''^"? will then 
represent the required subspace Xq. We call a set of values {Mfc}fc^,i the grid representative of 
the corresponding vector field. 

To obtain the element of X(r2) corresponding to m, we use a set of basis functions with ipk{x) 
representing a distribution centred around Xk. These basis functions are taken to have compact 
support and to satisfy the Partition- of- Unity (PoU) property 

Tig 

fc=i 
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The vector field X„ is then defined as follows: 

Definition 2.1. The vector field on Q whose grid representative is u takes the coordinate 
form 

k 

A plot of a typical basis function in one dimension is given in figure 01 



0.7 




Figure 3: Plot of a S-spline basis function i/'fc centred on ~ 7 with a grid width of 1. These 
basis functions satisfy the partition-of-unity property. 

Remark 2.2. In general, these vector fields do not commute amongst themselves in the Lie bracket, 
so they do not form a Lie subalgebra of X{^1). This will lead to a variational principle with 
nonholonomic constraints. Also in general, the value of Xu{xk) is not exactly equal to Uk, but is 
convergent to it in the continuum limit. 

Dynamics of a finite set of Lagrangian particles We shall proceed in describing our nu- 
merical method by introducing a finite set of Up Lagrangian fluid particles {Q/j}p'^i, whose 

velocities {Q/jj^^i are entirely determined by the grid velocity representation {uk}2'Li via the 
vector field as follows: 

Definition 2.3. The PoU vector field X^-.n^ G T^l"'' associated with a velocity grid represen- 
tative u is defined as 

X.;«,(Q)-EE"fc^'^(Q)-^' ge(Qi,...,Q„J=M'^^"''. (2) 

We shall constrain the dynamics of the particles so that a tangent vector Q may be represented 
as a PoU vector field evaluated at the point Q. That is, {Q, Q) lies in a distribution D^^^^ defined 
as follows: 



7 



Definition 2.4 (Tiie distribution D^^'^^). Let D^pm ^ j,^n^ distribution defined by 



I k 

Definition 2.5. ^ time series Q{t) (Qi(t), . . . ,Q„^(i)) wit/i G D^pm yto <t < ti 

is called a VPM trajectory. Each VPM trajectory defines a time series Uk{t) G M''^"9 x [io,^i] 
such that 

Q^(t) =^Mfe(t)V;fc(Q^), (3) 

fe 

for (3 = 1 , . . . , rip . This is the VPM tangent vector relation, which we will enforce as a 
constraint for the variational principle resulting in the VPM method. 

Remark 2.6. Given (Q, Q) one may invert equation 0) for the grid velocity representation u = 
. . . , ti„g) G K^'^^s X R, modulo the kernel oftpk{Qp) regarded as a matrix. (We shall see 
that this kernel does not affect the dynamics.) Later we shall write the Lagrangian as a function 
of Uk only and rely on this inversion to express the Euler- Lagrange equations for Q{t). We also 
note that a VPM trajectory Q{t) is specified entirely by {M/c(0}fc=i ^'^'^ initial condition Q{0). 
Changes of the initial conditions Q{0) for the VPM trajectories that leave invariant the grid velocity 
representation u{t) will provide the analog for VPM of "particle relabelling" in the continuum case. 



Gradient and divergence In this section wc describe how the operations of gradient, divergence 
and curl may be approximated using the particle-mesh discretisation. These approximations apply 
the two dual purposes of the basis functions ^pk- 

1. The tpk interpolate functions from the grid to the particles. 

2. The i/jfc also construct densities on the grid from weights stored on the particles. 

Notation: Square brackets [ • J*^ and [ • ]^ will denote these two maps from particles to grid and vice 
versa. Superscripts distinguish whether the quantity is evaluated on the grid or on the particles. 
That is, l']^ indicates mapping from grid to particles, and [•]'-' indicates mapping from particles 
to grid. 

Definition 2.7. Let {fk}^^! scalar quantity stored at the Eulerian grid points. Then 

k 

is an approximation of f evaluated at the particle locations. Furthermore, 

is an approximation of the gradient of the scalar f evaluated at the particle locations. 

Definition 2.8. Let {gpY^'^i be a distribution of values stored at particle locations. We construct 
a density on the Eulerian grid as 

I fS 

Furthermore, if the distribution is vector-valued g then 

I 13 

is an approximation to the divergence of g on the grid. 
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Discretised continuity equation Given a set of constant weights Dj3 on the particles /3 
1, . . . , np, to construct a density 

I 13 



one computes 



13 



and so 



I B k 



'^-[Df = -[V-{[ufD)f, 



so the corresponding grid representative [D]'~' satisfies a discretised continuity equation. 
Lagrangian for semi-discrete EPDiff 

Next we form the Lagrangian for semi-discrete EPDiff, as an approximation to the continuous 
EPDiff Lagrangian 

ic = ^^^(||M|P + a'||Vw|n dVol, 
in which the constant a has dimensions of length. 

Definition 2.9. Let {Nk{x)}2^j^ be chosen as a finite element basis so that functions may be 
approximated in the form 

/(a;) = ^iV,(a;)/fe, fk = fixk). 

k 

(This basis need not be the same as that used in the partition- of -unity representation of velocity.) 
Define the matrix H , which approximates applying the Helmholtz operator and integrating, as 

Hki = / Nkix)Ni{x) + a^VNkix) ■ VNi{x) d Vol, 

for some value of the constant a. Then the Lagrangian for discrete EPDiff is expressed in this 
basis as 

L{u) 2 X! ■ ^ki'^i = 2" ' ■ 

kd 

Remark 2.10. As in the continuous case, this Lagrangian is written entirely in terms of the 
Eulerian velocity (in this case, the velocity grid representation) . In the continuum case, this form 
of the Lagrangian admits Euler-Poincare reduction (as Eulerian velocity is invariant under the 
right-action of the diffeomorphism group Diff(ri)j. This reduction results in the EP equation 

Mt + ad* /X = , 

where fi = 6L/6u and ad* is the dual of the ad-action (Lie algebra commutator) of vector fields 
on the domain. In the VPM discretisation of EPDiff, an analogous equation will emerge, written 
on the Eulerian grid. 
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3 Variational principle for discrete EPDifF 



In this section we shall derive the equations for Q{t) from a variational principle applied to the 
Lagrangian and required to satisfy the VPM tangent vector constraint. Namely, the variational 
principle is constrained to restrict the solutions so that Q G (defined as the subspace 

{a : {a,Q) £ £)VPM| ^ Tq^I^^p). This constraint on VPM trajectories is the discrete analog of 
the Lin constraints in the Clebsch variational approach to continuum ideal fluid dynamics, as 
discussed for example in [HI ^2 • ^^le end of this section, we shall give a fully discrete variational 
principle which produces the numerical scheme. 



3.1 Constrained action principle for semi-discrete EPDiff 

We begin by defining the grid momentum as follows. 

Definition 3.1 (Grid momentum). In the same finite element basis as for definition \2.fA define 
the integration matrix M ( often called the mass matrix in the finite element literature ) as 



Mki = / Nk{x)Ni{x) dYol 
Jn 

The grid momentum m = (mi, . . . ,m„g) G ]R'^^"9 is then defined from the grid representative 
Lagrangian L{u) in via its derivative 

I 

This expression for the grid momentum is an approximation to 5L/5u in the continous case. 

Definition 3.2 (Constrained action). The action for semi-discrete EPDiff is defined in terms 
of three variables: the grid velocity u G IR''^"s; the particle positions G fi"''; and the Lagrange 
multipliers P p G T^ri"'' which will become the particle momenta on the Hamiltonian side. The 
action is given by 



A = i(w) + X^P/3- (^Q;ri-^WfcV'fe(Q/3)^ di. 



(6) 



This is the action for Lagrangian ^ when its particle velocities are required to satisfy the VPM 
tangent vector constraint given in 

Proposition 3.3. The variables (m, P, Q) which extremise the constrained action A in 0) satisfy 

= J2ukMQp), (7) 

k 

: Y^PPMQP)- (9) 



duk ^ 



Proof. After integration by parts, the first variation of A in [u, P, Q) is 



fi \ k / 
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and the result follows by direct calculation. □ 

Remark 3.4. [Left momentum map] Equation ^ in vrovosition UT^ bears a great resemblance 
to the momentum map for left action of the diffeomorphisms on embedded subspaces \1UI/ which 
describes the singular solutions of continuum EPDiff equation. We will see later that equation 0) 
is the discrete version of that momentum map. 

Remark 3.5. [Grid momentum] The grid-momentum relation \^ allows one to obtain u from 
dL/ dQ and Q by first calculating m, and then inverting the matrix H^i in 

dL 

This is the discrete analogue of the problem of solving for u from m in the elliptic relation 

SL , 
m = — = (1 — a" Ajit . 
ou 

for the continuous case jlOlj. Thus, the Lagrangian ^ is hyper-regular on the grid. 
3.2 Legendre transform 

We now pass to the Hamiltonian side via the Legendre transform, a process summarised in the 
following proposition. 

Proposition 3.6. The system of equations is canonically Hamiltonian with Hamiltonian 

function H given by 

^ ^ ^ E ^^-t'"^' ■ '"'^ ' (11) 

k,l 

where Uk = n H^^^Mmufn-a is defined in terms of Mkimi by inverting the matrix H^i in 
equation ^J^, and where rrik is obtained via equations ^ and 

Proof. We obtain the Hamiltonian via the Legendre transform 

H{p, Q) = E • E "^■^''^ (Q/j) - ^(") ' 

/3 k 

subject to equation Upon applying equation ^ the phase space action sum may be written 
as 

/3 k k (3 

= ^Uk-Mkivni, 

kd 

after switching the orders of summation. The Lagrangian Q may also be rewritten as 

L — —u ■ Au ~ —u ■ Mm . 
2 2 

Hence, proposition lH . 6l follows and we obtain the Hamiltonian 111(1 via the Legendre transform. □ 
Finally, we calculate Hamilton's canonical equations for this Hamiltonian. 
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Proposition 3.7 (Hamilton's canonical equations). Hamilton's canonical equations with H 
defined in equation j|) above may be expressed as 

k 

Proof. Hamilton's canonical equations are 



The P equation projects onto grid variables as 

. ^ _dH_^_^dH_ dmk 

Likewise, the Q equation projects onto grid variables via, 

• dH x - dH drrik 

" dP^" \-d^"dP^ 
= ^UkipkiQfj) , 

k 

as required. These are identical to equations (jSj) and ((T)), respectively. □ 
3.3 Constructing a fully discrete method 

To construct a fully discrete method we use the standard variational integrator approach as de- 
scribed in iI3J , applied to the constrained action principle in definition l3.2l We replace the integral 
over time by a Riemann sum over discrete time levels, and define the map 

: r?"" X T^T" 

which approximates Qp. We write the discrete action 

Aa^MY^i Liu-) - ^ m", Q""') 

n=0 \ P 

Minimisation of the discrete action over u, P and Q gives the numerical scheme. 
For example, consider the choice 



_ n-l 



yn /nri 1\ Q Q 



At 

In this case, the discrete action becomes 



N I 

Ad = AtY^r^HuiK-u^i 

n=0 \ kl 
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which is minimised by the solutions 

^ff,,«r = Y.p'^^'^iQVh (12) 

= Q^ + At^<+V.(Q^), (13) 

k 

= P»-AtP^+i.^««+i||i(Q^). (14) 

This system is equivalent to the 1st order symplectic Euler-A method (i.e. the 1st order symplectic 
method which is implicit in P and explicit in Q) applied to the Hamiltonian system given in 
proposition 13.61 

4 The discrete Euler-Poincare equation for VPM 

In this section we compute the discrete EPDiff equation directly on the Eulcrian grid. 

Theorem 4.1 (Discrete Euler-Poincare theorem). With the above notation and assumptions, 
let L{Q,Q) be a Lagrangian expressible as a function L{u{Q,Q)) of grid velocity representative u 
only. The following four statements are equivalent: 

(i) The VPM trajectory Q{t) is an extremal of the constrained action 



S= [ i(M)+^P/3- (Q^-^UfcVfc(Q/3)) dt 



/3 

with boundary conditions Q{a) ~ Q^, Q{b) ~ Q^,, and where Pp, P — 1, . . . , Up, are Lagrange 
multipliers. 

(ii) The VPM trajectory Q{t) is the solution to the canonical Hamiltonian system in proposition 
15'. 61 with suitable boundary conditions. 

(iii) The grid velocity u obtained from the grid momentum m using equation I^1U\) . itself obtained 
from the particle momentum P using lemma minimises the action S = ^^L{u) dt upon 
taking constrained variations 5uk which satisfy 

[H^ = <^'^k{t)ipk{Q 13) ^ '^k{Qfi)wk ~ [a.Au w]^, = [w\^ - [&du , 

k k 

where [ad^ii;]^ is defined by 

{ad^w)p = J2'^kiQp)(ukit) ■ ^^{Q^)^wi{t)~MQp)(wk{t) ■ ^^iQp))Mt) 



(iv) The grid momentum m ~ {mi, . . . ^mn^) G M''^""', satisfies the discrete Euler-Poincare 
equation 

rhfc + (ad* m)fc = 0, k = l,...,ng, 



whe 



(ad^m), = ^(A/-i),„(^^^(Q^)(P^.«,)^„(Q^) 



^J„i^P/3^^(Q^) • ^M,„(i)V'm(Q/3)), 
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so that 

(ad* m, w)g = (P, [adu w]^)p, 
where (•, ■)g is the grid inner product defined by 

k,l 

(i.e. a discrete approximation of the Lp' inner product in the continuous case), (•, •)p is the 
particle inner product on TQ^^p defined by 

{F, G)p = • G(3, 

and where Pp satisfies 

J2 Mkimi = ^ PpMQp)- (15) 

I 13 

Remark 4.2. The operation (aduK;)/3 in (iii) is the Lie bracket among vector fields evaluated at 
the particle location Q^. The operation a,d*^m is its dual with respect to the pairing (■ , ■)g on the 
grid. 

Proof, (i) <^ (ii) follows from the proposition 13. 61 

To prove (i) <=;• (iii) we note that the constrained variational principle given in proposition 13. 21 
is equivalent to the Lagrange-d'Alembert principle 



S dL dL 



6t dQg dQi, 



SQn = 0, 



with constrained variations SQ G i^vPM.Q and the constraint Q G £'vpm,q- The variations 6u 
must be expressed in terms of the variations 6Q and SQ which follows by taking variations in 
equation jSJ: 

dQ0 = .JQ^. (16) 

k k ^^1^ 

The variations 5Qp are written 

^Ql3 = ^Wk{t)ll)k{Q(3), 
k 

for some time series of velocity vectors on the grid Wk{t) which vanishes on the end points. Dif- 
ferentiating in time gives 

SQg ^Y.w,it)MQf,) + Y.^k{t)^{Qp) ■ Qp. (17) 
k k 

Combining equations p6|l and p7() gives 

^ 5uu{t)MQ,) - E MQp) Uk{t) + uk{t) ■ E ^{Qp)'^m ^k{t) ■ E W ) 

k k \ I I J 

which we denote as 



[5u]p = E Suk{t)iJk{Qfi) = E ^k(.Qp)wk - [adu w\p = [w]p ~ [ad^ w] 

k k 

This proves (i)<^(iii) and defines [ad^ii;]^. 
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Remark 4.3. The bracket-subscript notation introduced in the last formula emphasizes the 
VPM distinction between particle vector fields such as [it?]^ and their grid representatives Wk, 
related by [w]^ ~ Tlk '^k'^hiQ jj)- For example, 

(m, 5u)g = ^ rrikMki • = ^ ipi{Qp)Pi3 ■ Sui = ^Pi3 ■ [du]/3 = {P, [Su])p 

k,l 0,1 

where the momentum relation 0^ was used in the second step and the relation between VPM particle 
vector fields and their grid representatives was applied in the third step. A similar calculation allows 
one to write the dual relations defining the VPM particle- and grid-representatives of ad* . Namely, 

{P, [a.du w])p = (m, adu w)g , 
whose dual relation may be conveniently written as 

([ad* P]^, w}p = (ad* m, w)g , 
in order to define ad* in both particle and grid representations. In particular, this implies 

(ad* m, w)g = (P, [ad„ w]^)p , 

as claimed in the theorem. 

To prove (iii)4=>(iv) wc take variations Su in S: 

J (m, 5u)g dt 

\p, [Su])pdt 

{m,w)g - {P, [ad„ w])pdt 



= ^5 



{-fh,w)g ~ • (Mt)^iQp) ■Y.wiit)MQp) 

p \ I 



Wkit)^^{Qp)-J2Mt)MQp)] dt 



J {-rh, w)g - ^ Wk{t) ■ ^(g^)(P^ . ui)MQ0) 



1,13 



where we have integrated by parts. The grid representation it; is arbitrary and therefore 

+ ad* m 0, 
dt " 

as required. □ 
The correspondence (ii)<^(iv) was also proved by direct calculation in [S]. 
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5 Left action momentum map 



First, recall that a canonical action of a Lie algebra ^ on a symplectic manifold is a mapping 
from A to Hamiltonian vector fields on A4 which preserves the Lie brackets. Consider an element 
^ of yl and its action on which has Hamiltonian J. The momentum map J is related to 
the Hamiltonian J by 

for all such elements ^, where (■,■): A* x A ^ R is the inner product between A and and its 
dual A*. 

If A acts on a manifold A4 then we can define a canonical action of A on T*A4 with Hamiltonian 

This is called the cotangent lift of the action to T*A4. The definition of the momentum map for 
the cotangent lift of an action then becomes 

(J,0 = ((P,Q),Ia,)- (18) 
We define the left-action of X{^1) on il by 

d_ 
dQ' 

The Hamiltonian for the cotangent-lifted left action is then 

j(0(p,Q) = ((p,Q),^q). 



We wish to obtain a momentum map which maps into the representation of -Dvpm given by the 
map M^^"" X{n): 

k 

We do this by restricting J(^) to elements of -Dvpm: 

dQ 



J{u){P,Q) = l{P,Q)^Y,UkMQ)-^ 



\ k 
/3k 

and this relation defines the left action momentum map 

jf(P,Q)==^P^^fc(Q^). (19) 

As mentioned in Remark 13 . 41 this is again equation (jnj derived earlier from constrained variations 
of the VPM action ((HJ with respect to the grid representatives of the velocity. This momentum 
map is the discrete version for VPM of a general result for Clebsch variational principles for ideal 
fluid dynamics [51 lll|. 



6 Right action momentum map 

Right action of on {Q/3(i)}/3=i Next we define the right action of Diff(r2). To do this 

we require the entire "history" of u{t). Given initial conditions Q(0) = Q", the history of u(t) 
produces a solution Q{t) with 

Q(t) = ^Mfe(t)^fc(Q(t)), Q(0) = Q". 

k 
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This solution can be extended to a one-parameter family of diffeomorphisms gt '. Q with 

d 



dt 



9t{x) ^^Uk{t)ipk{gt{x)) , go{x)=x. (20) 



fe 



In particular, gtiQ'p) ~ Qpi^)- This allows a right action of ?/ e Diff(ri) on gt to be defined via 
composition: 

9t gt ■ V ^ 9t °V ■ 
Using the tangent map, one may define a right action of ^ e X(fi) on gt'. 



Again in particular, 



ox 



Differentiating equation (|20|l gives 

1 1 («) - E-c)^!^^ ^ t(«5) - S u.(.) t (OMO) ^ t («) . (2.) 

k k 

with initial conditions 

This means that, given Mfc(i) (and hence, given {Q/3(i)}^^^i), the Jacobian dgt/ dx{Q^p) may be 
obtained without needing to calculate gt as a map over the whole of f2. 

We can interpret this map as a canonical momentum map by extending the canonical coordi- 
nates {P/3, Qf3Yn=i to {-P/3; Q/3i '^/3}"=i where J/3 is a d X d matrix for each (3. We write J as a 
column vector, e.g., in two dimensions 

J = ('^ll,!: ^12,1, J21,l: J22,l, • ■ • , Jll,nj,, '^12, Hp, J21,np, ^22, Hp) , 

and consider the Hamiltonian system, cf. [5] 



where 



and 



Id 

-Id -B^{Q)K{Q,Jf 

K{Q,J)B{Q) 



BQ = u, if Q/3 = X] UkipkiQ/s) , 



dip 

ki 




When the Hamiltonian is a function of P and Q only, we recover the canonical Hamiltonian 
structure for P and Q. Furthermore, if the Hamiltonian is a function of grid momentum only, so 

that 

fe 



for some u. then 



as required. This larger system enables us to talk about discrete particle-relabelling, as summarised 
in the following theorem: 
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Theorem 6.1. Consider the time- continuous VPM discretisation of EPDiff, given in the above 
enlarged space, with Hamiltonian 



H = ^ MikTUk ■ {H ^)kiMijmj 



2 

klij 

with H the Helmholtz operator and M the mass matrix. Then the flows of the vector field with 
Hamiltonian 

13 

for any constant vector {^^j}^^^' ^s^'^s the Hamiltonian Ti invariant. 
Proof. The Hamiltonian h generates the flow 

Q = Jp-t 

P = -B^iQ)K^iQ,J)iP^), 

for P and Q, where we write 

(J • 1)^ = J/3 -C^, /? 1, . . . ,np, 

and 

= -P^/3^i./3' I3^l,...,np, i,j ^l,...,d. 
The Hamiltonian 7i is a fmiction of grid momentum m only so it suffices to check that the quantity 

^PptpkiQp), fc = l,...,ng, 
is invariant under the flow. We can check this directly, as 

= E (- {B^iQ)A^iQ,K)Pi)^MQp) + P(}^iQp) ■ Jp ■ ir)j , 

and the flrst term becomes 

~Y.^Jk{Qf3){B'^{Q)K^{Q,J)Pi)p = -Y.Klp{Q.J){Pi)0. 

13 (3 

showing that the momentum is invariant, as required. □ 
Corollary 6.2. The momentum map with 

•Jp = Pp ■ Jp^ 

is conserved for solution of semi-discrete EPDiff. 

Proof. The result follows directly from Noethcr's theorem, i.e., from invariance of the Hamiltonian 
TL in thcorcm l6.ll □ 

Remark 6.3. The symmetry which changes P and Q while leaving m invariant on the grid is our 
discrete form of the particle-relabelling symmetry. Next, we shall see that this symmetry results in 
a discrete version of Kelvin's circulation theorem. 
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7 Kelvin's circulation theorem for discrete EPDifF 

As discussed earlier, the discrete EPDifF Lagrangian is invariant under the right action of 
This means that in Corollary 16. 21 is a conserved momentum. In particular, 

=0, no sum on /3 

for each /3 = l,...,7ip. (This is obtained by integrating against a suitable function whose 
support contains only Qp.) 

We can interpret this result to prove a discrete form of Kelvin's circulation theorem. Consider 
a loop C{t) in which is embedded in the flow, i.e., 

C{t)^gtiC{0)). 



We choose C{t) so that some of the particles with trajectories Qpit) are located at the initial time 



< = on (5,3(0) G C{0). As gtiQ'3) = Qdit), those particles will stay on C{t) for all time. Define 



the set r so that (3 €T if Q° is located on C(0). 

In order to discuss the circulation theorem, we need to introduce a discretisation of density. As 
discussed in p], this is done by associating a constant Dp with each particle, so that the density 
on the grid may be written 

Dk=Y,DpMQp)- 

/3 

This allows us to represent m/D evaluated at the location of particle (3 as Pp/Dp. 
Next we need to approximate line integration round C{t). We do this by writing 



ic(t) D Jo Do-f^. ds 

where 74 : [0, 27r] C{0) is a parameterisation of the loop C{t). Substituting = 9t ° 7o yields 

Jc{t) D Jo Doj^ dx ds 

which we can approximate with a Riemann sum 

f m , \ - Pp ^ 

Jcit) D y Dp 

where 



Ax^ = Jp ■ -j—{sp)Asp, 

with 7o(s/3) = Qp, and Asp = sp+i - sp. 

Using this discretised line integration scheme, we can state our Kelvin circulation theorem as 
follows: 

Proposition 7.1. Let {uk{t)}2^-^ satisfy the discrete Euler-Poincare equations above, with {Dk{t)y^'Li 
satisfying the discrete density equation. Let C'(t) be a closed loop advected in the flow generated by 
the velocity 

u{x,t) =^Uk{t)il;p{Qp), 

k 

containing some subset of particles Qp, with /3 G B C (1, . . . Define the discrete circulation 



sum 



/36B Dp 
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Then I[t) satisfies 

Proof. The result proceeds directly from corollary 16.21 for the right action momentum map, which 
satisfies 

□ 

8 Numerical Results 

Convergence tests Wc begin by performing a convergence test for the ID equations 

mt + unix + 2mUx = 0, (1 — a^d'^)u, 

which, as discovered in is completely integrablc with the initial value problem dominated by 
peaked solitons (peakons) whose first derivatives are discontinuous. This property is illustrated in 
figure 21 which shows a numerical integration of the ID equations starting from a smooth initial 
condition with singular peaked solitons emerging in finite time. 




Figure 4: Numerical solution to the ID EPDifF equation with initial data u = (7r/2)e^ — 
2 sinhx arctan(e^) — 1, and scaling constant a = 1. The solutions arc obtained using the VPM 
method with 500 grid points, 1000 particles, cubic B-splincs for the basis functions, linear finite 
elements for the grid discretisation of the Lagrangian and a timestcp of 0.1 with the first order 
time discretisation given by equations p2N14(l . The figures show the velocity field at times 0, 5, 
25 and 50. At t = 5 the smooth initial condition has "leaned to the right" and a discontinuous 
peak has formed. By t = 25 the peakon is well separated from the smooth part of the solution and 
by i = 50 a second peak is starting to form. 

For our first convergence test we use the result given in [2] that for an initial condition u = 
(7r/2)e^ — 2 sinhx arctan(e^) — 1, with scaling constant a = 1, the asymptotic speeds of the emitted 
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peakons are 2/[(2n + l){2n + 3)], n ^ 0, 1, 2, . . .. In particular the asymptotic speed of the first 
peakon is 2/3. Figure |5l shows that the numerical calculation of the speed converges to the correct 
answer with a linear scaling for error against grid resolution. 



10" 




10" 



10' 

gridpoints per lengthscale 



10* 



Figure 5: Plot of error in calculating the asymptotic speed of the first emitted peakon from an 
initial condition u = (7r/2)e^ — 2 sinh x arctan(e^) — 1 against grid resolution (measured as number 
of gridpoints in one characteristic length a = 1), using the same method as figure 01 The number 
of particles used was twice the number of grid points, and the timestep used was scaled with the 
grid size to guarantee convergence of the fixed-point method for solving the linear system (i.e. 
not for accuracy). The logarithmic plot has slope -1 giving a linear scaling for the error with grid 
resolution. 

For our second convergence test we used the problem of an overtaking collision (illustrated in 
figure between two right-propagating peakons. gives a formula the phase-shifts for such a 
collision {i.e. the asymptotic difference in positions for the larger and smaller soliton with and 
without the collision). A plot of the error in the phascshift against grid-size is given in figure 13 

We found that the performance of the method when solving for head-on peakon/anti-peakon 
interactions was quite poor. During the collision the two peakons approach each other and stick 
together once they are both within a grid width of each other. This appears to be an issue with 
representing the momentum using Lagrangian particles, as to achieve a method with the correct 
results for the collision, the particle momenta would need to go to infinity during the collision. 
However, we can also view this as a benefit of the method. The peakon/anti-peakon solution 
represents an instability in the equations; namely that a small perturbation of the solution can 
result in a peakon/anti-peakon pair being created. As our method does not support this solution 
at the moment of collision, this type of instability does not pollute our numerical results. 

2D Flows 

In this section we show a few results obtained using the VPM method to discrctisc EPDiff in two 
dimensions with Lagrangian 



for a constant lengthscale a so that the velocity u is obtained from the momentum m = by 
inverting the modified Hclmholtz operator 




m = 



(1 - a^A)u. 
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10" 




Gridpoints per lengthscale 

Figure 7: Plot of error in calculating the fast shift in the faster peakon between a peakon of height 
1.0 and a peakon of height 0.5 against grid resolution (measured as number of gridpoints in one 
characteristic length a ~ I). This plot shows linear convergence of the error in the numerical 
solution. 
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In the first experiment the initial condition for the momentum liad a 2-dimensional "top-hat" 
profile 

1 ifa<x<6, c < y < d, 
otherwise. ' 



m = {m{x, y), 0), m 



for constants a,b,c, d, so that the velocity has continuous gradients and has compact support. 
Figure IHl shows the evolution of the velocity at subsequent times; it illustrates how EPDifF evolves 
to form singular filaments of momentum from smooth initial conditions. 




► 




Figure 8: Plots showing surfaces of velocity magnitude |m| at times t = 0, 0.45, 1.5 and 3.4 with 
a C^-smootli, compactly-supported initial condition for it. The momentum becomes supported 
on lines (so that the velocity has a "peaked" profile) which spread out. This is the 2-dimensional 
version of emerging peakons illustrated in figure^ These results were obtained using 65536 parti- 
cles, a 128 X 128 grid with periodic boundary conditions, tensor product B-spline basis functions 
and piecewise-linear finitelemcnt discretisation of the Lagrangian. The timcstep is 1.0 x 10"'^, and 
a = 0.2. 

Verification of conservation laws The next set of numerical results demonstrate the conserva- 
tion of the right-action momenta given in section and the connection with Kelvin's circulation 
theorem. 

Figure 1101 shows the value of the momentum map for right action Pp ■ Jp for a selection of 
particles from the flow in figure |H1 All the particles have the same value because they all have the 
same initial momentum and J is set to the identity initially. The figure shows that the numerical 
method preserves these conserved momenta up to round-off error. This follows from theorem 16. II 
for the time continuous equations and the fact that conserved momentum maps are also conserved 
by variational integrators. 
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Figure 9: Plots showing surfaces of velocity magnitude |m| at times t = 0, 1.55, 2.7 and 4.75 
showing an "overtaking" collision between two singular momentum filaments. The filament which 
is initially behind has greater momentum and so it catches up with the filament in front, transferring 
momentum to the front filament and causing a reconnection to occur. This is the 2-dimensional 
version of the process illustrated in figure This is the nonlinear reconnection process which is 
illustrated in the Space Shuttle image in figures d and 13 These results were obtained using the 
same method as figure |H1 
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Figure 10: Plots of Pp ■ Jp against time for a selection for particles with label illustrating that 
it is exactly conserved during the simulation. The upper plot is the x-component and the lower 
plot is the y-component. Any variation seen is due to numerical round-off error. 
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To verify the discrete circulation conservation discussed in sectional we took an arbitrary loop 
containing some of the particles and advected the loop with the flow shown in figure |21 using 



Q"+^(,s) = Q"(,s) 



k 



UkipkiQis)), 



where s parameterises the loop. During the course of the flow this arbitrary circulation loop 
evolves, changing shape and length significantly. However, the circulation around the loop remains 
constant (up to numerical round off), as verified numerically in the following. 






Figure 11: Plots showing an embedded loop in the time- varying flow obtained from a solution of 
EPDiff, illustrated in figure|Hl at times (left-to-right, top-to-bottom) t = 0, 0.25, 0.5, 0.75, 1, 1.25, 
1.5, 1.75, and 2.0. 

To write down the circulation integral, we choose an initial density 

Dk = {M-^)kiY,DpMQp). 

where the values of D do not matter much, as they are not coupled with the dynamics. Hence, we 
choose the values Dp = 1, /? = 1, . . . , Np. To obtain the discretised loop integral 

we need to calculate Aa;/3 as discussed in section [7| This is done by finding discrete line elements 
Aa;° for the initial loop and then calculating Acc for subsequent timesteps using 



Summary of circulation loop figures 
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• Figure ITTI shows the evolution in time of this circulation loop, under the flow induced by the 
expanding waves in figure |S1 

• A plot showing initial and advected line elements is given in figure E| This plot illustrates 
how the line elements evolve when a loop is stretched out by the flow. In the top and the 
bottom of the loop, where the stretching is greatest, one can see how the line elements extend 
to provide a numerical approximation to dx on the loop. 

• Finally a plot of the circulation integral is given in flgure This plot shows that the 
circulation round the loop is exactly preserved during the simulation (up to round-off error 
in the calculation of the discrete integral). 

9 Summary and Outlook 

In this paper we studied the Variational Particle-Mesh method applied to the EPDiff equation. 
We introduced a constrained variational principle for the method and gave discrete Euler-Poincare 
formulae on the Eulerian grid resulting from the variational principle which show that the grid 
velocities and momenta satisfy the EPDiff equations in Eulerian form. Next we looked at left- 
and right-actions of velocity vector fields on the Lagrangian particles and obtained corresponding 
momentum maps. The left-action, when restricted to the finite space of velocity fields used in the 
method, gives rise to a momentum map which is the same formula as used for calculating the grid 
momentum from the particle variables. The right-action, which had to be interpreted in a wider 
space, can be interpreted as a discrete form of particle-relabelling since it corresponds to moving 
the particles and also changing the momenta in such a way so that the grid velocities remain 
constant. Finally we gave some interpretation of these transformations in terms of matrices which 
determine the local deformation of infinitesimal line elements, thereby allowing us to write down 
discrete loop integrals on advected loops. This led to a discrete circulation theorem. 

Our next aim is to find an extension of this work which gives a discrete circulation theorem 
for fluid PDEs which involve mass density and other advected quantities as well as velocity. The 
general approach, following the continuous theory, will be to 

• specify the transformations corresponding to discrete relabelling, 

• determine the transformation of density and other advected quantities under this discrete 
relabeling group, 

• calculate the momentum densities obtained from these transformations. 

• show that the ratio of these momentum densities to the mass density is invariant. 

Including advected quantities in this awy will allow introduction of potential energy and hence 
linear dispersion effects into the numerical description of the internal wave interactions using the 
VPM method. 
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Figure 12: Plot showing the curve with embedded line elements at time t = and t = 2.3. The 
line elements at time t = 2.3 are the exact tangents to the curve which passes through the tangents 
given at time t = and is then advected using the time independent flow given in figure |H1 
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